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Abstract 

As known, any numerical simulation is composed of two parts: (1) the initial part of 
writing the relevant code and (2) the running of this code on the computer screen. The 
second part of running the program is extensively discussed theoretically and technically 
in the relevant literature. In this work we pay special attention to the less discussed 
first part and show that it may be discussed in a terminology and notation which de- 
scribe physical phenomena. As examples we discuss the two cases of simulating (1) the 
harmonic oscillator and (2) the electron- photon interaction which results in the known 
Lamb shift. 
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1 Introduction 

The problem of validating scientific theories through numerical simulations have been dis- 
cussed from several points of view [H, EJ [3]. There is now almost no scientist (physicist, 
chemist, biologist etc) that does not use the powerful means of numerical simulations as a 
necessary tool in his research. Moreover, it is accepted [2, [3j that if some scientific theory 
is found in its numerical simulation version to be valid on the screen then generally it is 
valid also outside it. Thus, a correspondence may be drawn between the various stages 
in physical theories of first proposing the scientific theory (writing the relevant equations), 
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and then testing it through experiments to the analogous steps in computer simulations of 
first writing the program and then running it on the computer screen. This correspondence 
between the two processes is especially emphasized in the experimentation stage except for 
the differences due to their different characters. That is, whereas the physical theories are 
proved or refuted through real experiments performed in the (three-dimensional) laboratory, 
the relevant "experiments" in the numerical simulations are the running of the involved 
programs on the (two-dimensional) computer screen. 

First of all we note that a programmer who wants to numerically simulate any process 
may accomplish his task by writing any of a large number of different codes which are all 
equivalent for obtaining the same simulation. Also, when he begin to write his program there 
is no way to predict beforehand the final code (from all the possible ones) so the writing 
stage of any numerical simulation has some unpredictability related to it. 

In this work we are, especially, interested in the programming procedure and in the 
possible mathematical description of it. For this purpose we use the known fact that the 
written code of the numerical simulation of any process must contain the full and detailed 
description of the simulated phenomena. That is, only when one introduces into the code in 
advance all the details of the simulated system that he may expect an appropriate simulation 
of it when this program is later run on the computer screen. Thus, since the programmer 
must be very careful in correctly describing in his code the behaviour of the simulated 
system one may suppose that the code-writing itself evolutes in such a manner that it 
reflects the real evolution of the simulated process. That is, referring to the correspondence 
between numerical simulation and physical theories, we assume that as the latter describe 
real processes as rates of change of some variables in the infinitesimal (or infinite) limits of 
other dependent variables so the simulation of this process may be described in an analogous 
way. 

One may analyze the numerical simulations from the point of view of two interacting 
systems; the simulated system and the one that enables, through writing the relevant code, 
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the numerical simulation of the former. The dominant and active system in this supposed 
interaction is of course the simulated one which actually dictate the content of the written 
code which should, as remarked, represent and describe the simulated system. Thus, using 
the terminology of Synergetics (see, for example, page 195 in [4]), the programming act is 
"enslaved" by the simulated system and may be actually discussed in terms of it. This is the 
basic and central property of Synergetics which have been successfully applied not only to 
Physics bu also to Chemistry, Biology and other exact disciplines. That is, one may describe 
the interaction between any two heavily interacting systems from the exclusive standpoint of 
the dominant system. In other words, the evolution and development of the "enslaved" [4] 
and passive process of writing the code may be entirely discussed in terms of the "slaving" 
and dominant simulated system as done here. 

Thus, for taking into account the fact that there are many possible different codes which 
fulfill the same numerical tasks we introduce, as done in other analogous situations which 
involve different possible ways to obtain the same result (see, for example, [HIE]), an extra 
variable. We use here the known Parisi-Wu-Namiki Stochastic Quantization (SQ) method [TJ 
8J by which an additional variable have been introduced into either the Langevin equation [9] 
or the Fokker-Plank one |10| . In this formalism one assumes that some (generally unknown) 
stochastic process [1 1 J occurs in the extra dimension of the additional variable and the 
equilibrium physical situations are approached in the limit of the elimination of this variable 
which is obtained by equating all its different values to each other and taking to infinity 
[HIE]- This formalism is appropriate for the program- writing stage of the simulation which 
may evolute along different equivalent routes and, therefore, the unpredictable element of it 
corresponds to the mentioned stochastic process. The numerical equilibrium stage is obtained 
when the writing process ends and one remains with the finished program. 

For a suitable analysis of the possibilities allowed at the code-writing stage we discuss 
a large ensemble of programmers which all try to perform the same numerical simulation. 
We calculte the probability to find all or most of them obtaining not only the same final 
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simulation but also writing the same code which truly represents the simulated phenomena. 
That is, we want to find the path integral [T2] correlation and the conditions which must 
be fulfilled in order to obtain a large probability for finding them with the same numerical 
evolution which really describes the simulated system. Thus, our ultimate test will be to 
find out if at the equilibrium stage after they have finished writing their codes the obtained 
simulation correctly represents the simulated real systems. This indeed will be shown in the 
following sections by the obtained expression for the correlation between the ensemble of 
codes which turns out to be very similar to the correspnding expression for an ensemble of 
the real simulated systems. 

In Section 2 we present the SQ method and relate it to numerical simulation. In Section 3 
we restrict the discussion to the computer simulations which simulate the physical harmonic 
oscillator. We show, by discussing, as remarked, the simulation process in terms of the 
simulated harmonic oscillator, that the expression obtained for the correlation between the 
computers of the ensemble are very similar to the known expression for the correlation 
between the members of an ensemble of real quantum harmonic oscillators. 

In Section 4 we discuss the numerical simulation of the known Lamb shift process [EH 
[HI [15] and show, using the SQ formalism, the Fokker-Plank equation [10], and the classical 
Feynman diagram [El [EH \T7\ that at the numerical equilibrium state one may obtain for the 
correlation, as for the harmonic oscillator case, the analogous known expression obtained 
in the framework of quantum field theory [T4l [18]. In Section 5 we further analyze and 
summarize the discussion. 

2 Representation of the simulation process as a stochas- 
tic Langevin equation 

As remarked, the code of any simulation process must be written so as to include and describe 
all the details of the simulated system if one wants the simulation to truly represents it. 
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Thus, one may expect that the programming act evolutes and developes in a manner which 
reflects the corresponding evolution of the simulated systems. And since, as noted, these real 
systems are mathematically described by taking rates of changes so the numerical simulation 
of them may likewise be discussed in such a manner. One have only to take into account 
the additional functionality, taken care of by the extra variable (denoted s) , that by running 
different codes one may obtain the same numerical result. Thus, if we denote the process of 
writing the code by the programmer i by where 1 < % < N, then these qj may be analyzed 
by taking the rates of their changes with respect to s, that is, as the generalized Langevin 
equation [9] 

= Ki(q( 8) t, x)) + Vi(s, t,z), i = 1, 2, . . . N, (1) 

where N denotes the remarked A-member ensemble of programmers which all try to perform 
the same numerical simulation. The t^'s denote stochatic processes in the variable s. As 
remarked, these processes stand here for the described unpredictable nature of the program- 
ming act where one can not predict beforehand the final code. The variables % depends 
upon s and upon the spatial-time axes (x,t), where x denotes the two dimensional spatial 
axes of the screen since these numerical simulations occur on the screen and t is the time. 
The Ki are given in the SQ theory by [3 [8] 

dS- \q] 

Ki(q(s,t,x) = -(—^-) q=q ( s ,t,x), (2) 

where Si are the actions Si = J dsLi(q,q) that determine the forms of g» and Li are their 
Lagrangians. In order to discuss the "evolution" of the code- writing process which determines 
the simulation on the screen, we consider the time and s intervals (to,t), (s , s) and divide 
each of them into N subintervals (to, ti), (ti, t 2 ), ... (£jv-i, t) and (sq, si), (si, t 2 ), . . . (sn-i, s). 
The subdivision of paths into small subintervals is basic and central in the path integral 
method not only in the quantum version of it but also in its application to classical systems. 
This subdivision is very important here for the two discussed examples of the harmonic 
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oscillator in Section 2 and the Lamb shift in Section 4. Since only by applying and calculating 
the relevant expressions (the Langevin equation here and the Fokker-Plank equation in the 
next section) in these subintervals the results shown in the following are obtained. We assume 
that the Langevin Eq (DQ) is satisfied for each member of the ensemble at each subinterval 
with the following Gaussian constraints [8] 



where the angular brackets denote an ensemble average with the Gaussian distribution. The 
r signifies the N subintervals of each member and the i, j denote these members where 
N > i j > 1. Since, as remarked, the written code should fully describe the evolution of 
the specific simulated process the a have different meanings which depend upon the identity 
of this process and the context in which Eqs (PQ) and ([3j) are used. Thus, in the classical regime 
a is [8] a = ^y-, where kp, T, and / are respectively the Boltzman constant, the temperature 
in Kelvin units and the relevant friction force. In the quantum regime a is identified |8j 
with the Plank constant h. We note that using Eq (pQ) together with the specific constraints 
from Eq j3]) enables one [8j to discuss a large number of different classical and quantum 
phenomena. It has been shown [8] that the right hand side of Eq ((3l) may be written as 



<r)i(t r ,s r )>= 0, <r)i(t r ,s r 



)rjj(t r , s r ) >= 2a5ijS(t r — t r )S(s r — s r ), 



(3) 




1 



exp( 




(4) 



which is the probability to have a value of r\i in (yi,yi + dy) |8j, where 



Vi = 



dqj(s,t,x) 
ds 



Ki(qi(s,t,x)) 



(5) 



We want to calculate the Green's functions (t 1? t 2 , • • •) which determine the correlation 
among the members of the ensemble [8j. Thus, one may define, as in j8], the Green's functions 
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which depend also upon the variable s. 



\j...{t Q , s ,h, sx, . . .) =<qi(t , s )q j (t 1 , Si) . . .>= (6) 

Si(q(t,s)) 



C J Dq(t, s)qi(t , s )q j (t 1 , Si) . . .exp( 



a 



where Si are the actions Si = J dsLi(q,q), C is a normalization constant, and Dq(t,s) = 
Ill=i dqiit, s). As seen from the last equation the Aij (ti,t 2 , ■ ■ ■) were expressed as path 

iS(q) 

integrals [12]. Note that the quantum Feynman measure e n is replaced in Eq (p) and in 
the following Eq j7]) by e~~^ as required for the classical path integrals [H [IU [20] . 

It can be seen that when the s's are different for the members of the ensemble so that 
each have its specific Si(q(si,t)), Ki(q(si,t)), and r]i(si,t) the correlation in j6]) is obviously 
zero. Thus, in order to have a nonzero value for the probability to find a large part of the 
ensemble writing the same numerical code we have to consider the stationary configuration 
where, as remarked, all the s values are equated to each other and eliminated. For that 
matter we take account of the fact that the dependence upon s and t is through q so this 
ensures [8] that this dependence is expressed through the s and t differences. For example, 
referring to the members i and j the correlation between them is Ay(^ — tj, Sj — Sj), so that 
for eliminating the s variable from the correlation function one equates all these different s's 
to each other. We, thus, obtain the following stationary equilibrium correlation [8] 

A ij „Xto,s ,...) s t=<qi(t ,s )q j (ti,s 1 ) ...> st = C J Dq(t)qi(t )qj(ti) . . . (7) 



exp(- 



where the subscript of st denotes the stationary configuration. In other words, the equilib- 
rium correlation in our case is obtained when all the different s values which give rise to 
different programs and so to different simulations are equated to each other in which case 
one remains with the set of programs which simulate the same system. 
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Thus, if all the members of the ensemble write the same version of the program which is 
mathematically described as having similar actions S (in which the s values, which denote 
here different codes, are equated to each other) one finds with a large probability these 
members, in the later equilibrium stage, with the same code. That is, introducing the 
same similar actions into the corresponding path integrals one finds this mentioned large 
probability. This is exemplified in the following Section 3 for the harmonic oscillator and in 
Section 4 for the energy shift case. 

3 The numerical simulations of the harmonic oscillator 

As an application of the former discussion we calculate the correlation of the N member 
ensemble with respect to the numerical simulations of the harmonic oscillator. That is, we 
calculate the probability to find all or most of the programmers writing the same code of 
the harmonic oscillator which truly represents the harmonic oscillator. At the initial time 
t when all the programmers just begin to write their respective versions of the harmonic 
oscillator code they are certainly identical to each other. That is, all the programmers begin 
from the same common starting point (denoted q°) and we want to find the probability 
that they end at the later time t with the same code (denoted g( 2Ar-1 )) which refers to the 
harmonic oscillator. We assume, for convenient mathematical representation of the following 
discussion, a 2N member ensemble and also a subdivision of each of the intervals (t , t) and 
(s , s) into 2N subintervals (t , ti), (tt,t 2 ), • • • (W-i^) and (s , si), (si, s 2 ), . . . (s 2 n~i, s). We 
write the Langevin equation (Op) for the subintervals (tk-i,tk) and (s^-i, Sk) in the form [8] 



q\ - qt 1 - K t {q k - l ){s k - 8k . x ) = drf" 



(8) 
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The appropriate Ki for the harmonic oscillator is [8] 



Kitf-Ht k , a k )) = m|(^) - (^p-) q=q(t , s) , (9) 



where ~ ^ _ t » . The c%(s) are conditioned as [8] 



<dr)i(s) >= 0, <dr]i(s)dr]j(s)>- 



for s 7^ s 



2aSijds for s = s 



where the a is as discussed after Eq (|3j) and the probability from Eq (jlj) assumes the following 
form for the harmonic oscillator [8] 



P(g fc , 4, s k \q k ~\ t k . u s k -!) = ( , ) 2N ■ (10) 

'27r(2a(s fc - s fc _i)) 



Pl Y 2(2^,-^0) J ' 

which is the probability that the drfc _1 from the right hand side of Eq (jSJ) take the values at 
its left hand side |8j and the index % runs over the IN members of the ensemble. Note that 
Eqs fl8l)-(fT0l) are discussed in [H] without relating the variable s to any numerical simulation. 
Here, we relate the variable s to the possible numerical codes of the harmonic oscillator 
which are written, as remarked, so as to reflect and describe all the details of it. Thus, the 
evolution of the code- writing process should reflect that of the simulated harmonic oscillator 
according to Eqs (l8l- (fT0l) . This is the meaning of saying that the right hand side of Eq (jHJ), 
which represents the unpredictability of the code-writing process, should truly reflects the 
left hand side of it which represents the harmonic oscillator. A Markov process [TTJ in which 
rj(s) does not correlate with its history is always assumed for these correlations. Eq (TTOl) is 
the probability that the ensemble is found at t k and s k with the harmonic oscillator code 
q k if at tfc_i and s k -i it was at the similar q k ~ 1 . Note that by the word code we do not 
necessarily mean the complete program of the harmonic oscillator, that is, we use this word 
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even for a very small part of it. The probability for the entire interval that the ensemble is 
found at t and s to be with the same harmonic oscillator code g( 2iV_1 ) if at the initial to, so 
it was at q° is [8] 

P^" 1 ), t N , s N \q°, to,3 ) = J J J P(g( 2JV - 1 ) ) t N , s N \q^ N - 2 \t {2N . 2) ,s {2N ^) ■ ■ ■ 
■ ■ ■ P(q\ t k , s k \q k -\ t k -x, s fc _i) ■ • • P(q\h, Sl \q°, t , s )dq {2N ^ ■ ■ ■ dq k ■ ■ ■ dq 1 (11) 

In order to be able to solve the integrals in the former equation we first substitute from Eq 
(111) into Eq ((SJ). Thus, dividing the result by the infinitesimal interval s k — s k -\ = 5s, writing 
for V(q) the quantum mechanical potential energy V(g») = \mwoq 2 with the eigenvalues 
E v — Wq{v + =) v — 0, 1, 2, . . . and Fourier transforming we obtain for Eq (jHJ) [8J 



dqj(n,s) 
ds 



-m((K k ) 2 + (w ) 2 )qi(Kk, s k ) + fji(Kk, s k ), (12) 



with the following Gaussian constraints (in which we denote the Fourier transforms of qi(t, s k ) 
and r]i(t, s k ) by s k ) and fji(K k , s k ) respectively). 

<Vi( K k,s k )>= 0, <fj i (K k ,s k )f} j (kk,s k )>= 25 ij 5(K fc + k k )5(s k - s k ) 

In the following we write s k -i for s k - Solving Eq f fl2l for qi(n k , s k ) one obtains [8] 

qi(K k , s k ) = q exp(-m((n k ) 2 + (w ) 2 )s k ) + / exp(-m((/? fc ) 2 + (13) 

Jo 

+ (w ) 2 )(s k - S k -l))fji(K k , Sfe_i)rfs fc _i 

From the last equation we obtain the correlation Py 5 («fc, s k — s k -i) 

Pij s {K k , s k - sjfc-i) =<q'i(K k ,s k )q'-(K k ,s k --i)>= (14) 
exp(-m((K fc ) 2 + (wo) 2 )|sfc - s fc _i|) 



M( K k) 2 + (u>o) 2 
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Since we want our results to include a time dependence we Fourier transform Eq (PHI) back 
to obtain 

If 1 

PijAtk -t k -i,s k - Sfc_i) = — / dK k —- — - T — — —exp(iK k (t k - t k -i) - (15) 

H Air J ml I Ki. r -+- l?/v>l z l 



2it J m((K k ) 2 + (w ) 
-m{{K k ) 2 + (w ) 2 )\s k - sjfc-il) 



The former equations (fl2j) - ([T5j) are for the subintervals (t k _i,t k ) and (sfe-i,Sfc) so that for 
obtaining the corresponding expression for the whole intervals (to,t) and (s ,s) we use the 
following property of correlation functions [19] that if <Qi(xi)qj(xj) >=<^(xj) ><qj(xj) > 
then 

k=N-l 

<qi{x 1 )q 2 (x 2 ) . . .q 2 N{x2N>= II <( ?2fe+l(2:2fc+l)92fc+2(^2fc+2) > 

fc=0 

Thus, the generalization of Eq ( TT5l ) to the entire intervals is 



P{q {2N 1) ,t,s\q°,t ,s ) st =<q (t ,s )qi(t 1 ,s 1 ) . . . q 2N -i{hN-u «2J\r-l) >= 
j fe=JV-l j 

= A 7 ^ m( , )2 + ( , 2) exp(iK fc (t 2fc+1 - t 2k ) - (16) 



fc=0 
\2 1 /„,. \2 



-m((« fc ) + (lW ) )|«2fe+l - «2fc|) 

Eq (fTBl is the sought for probability to find at t and s the whole of the ensemble of 2N 
programmers having the same harmonic oscillator code g 2 ^" 1 if at the initial to and So they 
all begin from q° which is the initial stage at which they begin to write their code. As 
remarked, the stationary configuration is obtained in the limit of eliminating s so equating 
all its different values to each other, as required by the SQ method, one have 



k=N-l f ,-w \t 2 k+i-t2k\ 

n 



k=0 



k=0 
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Note that the elimination of the variable s is obtained by only equating all its different 
values to each other without having to take the infinity limit. Figure 1 shows the correlation 
from Eq (fTTt l as a function of t for m = 1 and wo = 0.4. It begins from the value of 1.25, 
which corresponds to the assigned values of m and wq, then steps through a maximum and 
vanishes for large t. Figure 2 shows a three-dimensional graph of the general correlation 
from Eq ( JT6l ) as function of t and s and for the same values of m — 1 and wo = 0.4 as in 
Figure 1. Note that for large s the correlation vanishes even at those values of t at which it 
attains its maximum in the stationary case of Figure 1. The values of m = 1 and wo = OA 
are typical values used for simulation purposes. 

Note that the stationary state from Eq (fTTI ) have been obtained by inserting the har- 
monic oscillator Langevin expression from Eq (fl2l) into the action S of each subinterval pair 
(tfc_i,tfc), (sfc_i,Sfc), 1 < k < 2N of each member of the ensemble. This kind of substi- 
tution is clearly seen in Eq (fTOl which includes the Langevin relation from (jHJ) in each pair 
of subintervals (tfc-i,tfc), (sk-i,Sk). Note that the substituted expressions of the harmonic 
oscillator into the actions of the subintervals are, of course, not identical since in this case 
the probability to write the same code by all the programmers is trivially unity. As one 
may realize from Eq (fT6l) the substituted expressions differ by s and t and only in the limit 
that these expressions have the same s and t that one finds the same code shared by all the 
ensemble members. 

As noted, the ultimate test of any simulation is that it correctly describe the simulated 
phenomena when the written code is run on the computer screen. Thus, taking into account 
that the classical path integrals discussed here are formulated in the Euclidean formalism 
[TBI 120] in which the time t is imaginary we see that the correlation from Eq (fT7l) is almost 
the same as that of the quantum harmonic oscillator which is [8] 



k=N-l p -iw \t 2 k+i-t2k\ 



A, 




t2N-2) — Yi 



<iiimw 



(18) 
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Thus, as remarked, introducing into the numerical code of the harmonic oscillator all the 
detailed description of the real one we obtain the correct simulation as seen from comparing 
Eq f TTTT l to (1T8|) . Moreover, when the remarked substitutions of the harmonic oscillator 
relation is performed in a dense manner over very short intervals of t and s, in which case 
the substituted expressions are almost identical, one may obtain the situation in which all 
the members of the ensemble write exactly the same code and, therefore, the correlation f fT7l ) 
becomes large. In such case iV becomes very large and may be written as N = where 
St is the time difference of each of the 2N subintervals. Thus, Eq ( fT7l ) becomes 

p —NwoSt p —woSt 

^-S M | A ^^ = I5 ^ = (_)» (19) 
From the last equation one realizes that if the condition 

e -w 8t = 47Tm w 0} (20) 

is satisfied then the correlation among the ensemble members is maximal because each of 
them writes exactly the same code so the mentioned probability is unity. Note that in this 
case not only the s intervals tends to zero but also the fs as seen from the former equations. 
In this case the left hand side of Eq (1201) becomes almost unity (for not very large values of 
Wq) and thus for having a probability of unity one have Wq = t~ (keeping the former value 
of m = 1). 



4 The Lamb shift example 

We see from the former section that substituting the harmonic oscillator expression into the 
actions S of the path integrals [T2], [16] which are related to the large ensemble of programmers 
establishes it among them in the sense that the probability to find them writing the same 
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program which represents the harmonic oscillator is large. We show this again for the 
example of a two-state electron which emits a photon and then reabsorbs it where the total 
energy is not conserved. We assume, as for the harmonic oscillator example, that there are 
a large number of different numerical codes of this process which reflect the large number of 
different ways which lead to the same simulation. Thus, as for the harmonic oscillator case, 
we introduce a large ensemble of programmers all of them want to simulate this electron- 
photon interaction upon their computer screens. As in the former section the differences 
among the written programs may be related to the different values of the extra variable 
s. And, as for the harmonic oscillator case, the equilibrium state is obtained when all the 
different values of s are equated to each other and taken to infinity. 

We subdivide the intervals (sq, s) and (to, t), during which this process occurs, into a large 
number of subintervals (sq, si), (s±, S2), ... (sjv-i, sjv) and (to, £1), (£1,^2), • • • , (£jv-1j ^n) an d 
formulate the appropriate expression for the described electron-photon interaction over the 
representative subintervals (t k _i,t k ) and (s k ~i,s k ). We calculate the probability to find the 
ensemble of programmers writing the same code which truly simulates the remarked electron- 
photon interaction. We find it better to discuss now this probability using the Fokker-Plank 
equation [HI [10]. That is, we begin from [8l HP] 



where P(q k , t k , Sfe|9 fe-1 , £fc-i, Sfc— 1) is given by Eq (TT0~|) and denotes the probability to find the 
relevant ensemble of programmers at t k and s k having the code q k if at the former t^-i and 
they have the code q k ~ x . Note that, as in the harmonic oscillator example, the word 
code does not necessarily means the complete program of this process, that is, even a very 
small part of it is called code. The operator F(q k ) is [H] 



dP(g k ,t k ,s k \g k 1 ,t fc _i,s fc _i) 
ds 



F(q k )P(q k , t k , s k \q k \ t k - U s fc _ x ) 



(21) 



F(q k ) = H(q k ,7r k ) 



(22) 



THE LAMB SHIFT EXAMPLE 



15 



where H and n k are the "stochastic" Hamiltonian and momentum respectively and a is as 
discussed after Eq j3]). The momentum operator 7r fc is defined as in quantum mechanics [8 J 
■n k = —2a-£pr, and its commutation with the operator q % satisfy [8J [7r fc ,(f] = 2aSti, where 
all one have to do in order to obtain the quantum regime is to set [8] a = %. From the 
former relations one may develop, as has been done in [8], an operator formalism similar to 
that of quantum mechanics, especially, the corresponding "Schroedinger", "Heisenberg" and 
"interaction" pictures. 

Using the former discussion we find the conditional probability to find at s and t the 
ensemble having the code q N if at the initial s and t they have the code q°. That is, 
one can write this probability in the "interaction" picture for the intervals (t, to), and (s, So) 

mm 

P\q N , t N , s N \q°, to, so) = P 1 \q° , t , s ) + (23) 
+ J F I (q N )P I (q N - 1 ,t N . 1 , SN . 1 \q°,to,So)dq N - 1 , 

where P 1 (q°, t , So) is the probability to find the system at the initial to and so with the code 
q° which probably contains at the initial stage only a few bits and the superscript / denotes 
that we consider the "interaction" picture. Note that q depends upon s and t so the integral 
with respect to q is, actually, a double one over s and t. Substituting, in a perturbative 
manner [12j for P\q N ~ l , tjv-ij sn-i\q°, to, Sq) on the right hand side of Eq (l23l one obtains 

n=oo -. q N N N 

pi(q",t N ,s N \q°,to,So)=j:^ L L ■ ■ ■ /, dq N T (F 1 ^ 1 ) F 1 {q 2 ) . . . 

n= Til Jq° Jq° Jq° 

...F I (q N )P I (q ,to,So) = P I (q°,to,So)+ I dq 1 F 1 (q 1 )P I (q° , t , s ) + (24) 

J q° 

+ dq 2 f dq'F\q')F\q 2 )P\q\ t , s ) + . . . f dq N ~ l f ^ dq N ~ 2 . . . 

Jq° JqO Jq° Jq° 

... f q F I (q 1 )F I ( q 2 )...F I ( q N )P I (q ,to,So) 
We, now, follow the same rules in [13], except for the introduction of the variable s, for 
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representing the electron and photon before and after their interaction. The extra variable 
s is introduced into the relevant quantities so that in the limit of s — » oo, as required in 
the SQ method [HE], the known expressions [13] which represent the electron and photon 
are obtained. We note first that the probability P 1 is no other than the state of the system 
[8] (as in quantum mechanics the system states of the SQ theory have also a probabilistic 
character). Thus, in the former equations we may assign to the initial So and to the value of 
zero and refer to -P 7 (g°, t = 0, s = 0) as the initial state of the ensemble system. This initial 
common state denotes, as for the harmonic oscillator case discussed in the former section, the 
common starting point of all the ensemble programmers which simulate the electron-photon 
interaction. 

As remarked, the electron is assumed to have two different states so that at t\ and s\ 
it was at the higher state 2 from which it descends to the lower one 1 through emitting a 
photon. Then at t 2 and s 2 it reabsorbs the photon and returns to state 2 as schematically 
shown at the left hand side of Figure 3. The incoming electron and the emitted photon at 
ti and Si may be represented by e~ le2tl + e - 4e 2Si(i-i<5) an j e -tw x ti _|_ e -vw x si(i-tS) respectively, 
where 5 is an infinitesimal satisfying 5 ■ oo = oo, and S ■ c = 0, (c is a constant) pEE]. This 
is done so that in the equilibrium configuration, which is obtained in the SQ theory when 
s — > oo, the terms in s vanish as required [H [8] and one remains only with those in t as in 
[IB] . The outgoing electron after emission at t\ and s\ may be represented by the plane wave 
e ieiti _|_ e ie lSl (i+i&) wriere § has the same meaning as before. At the reabsorption stage at 
t 2 and s 2 the electron is represented, before absorbing the photon, by e~ leit2 + e - le ^0-~ l5 ) 
and after the absorption by e l€2t2 + e * €2S2 ( 1+i<5 ). The photon is represented at the reabsorption 
stage by e lWxt ' 2 -\- e tWxS2 ^ 1+tS \ Also, the emission itself, denoted by the vertex in Figure 3, may 
be represented, as in the quantum analog [H], by g\ a and the reabsorption by g^ s , where 
an explicit expressions for g\ s and g^ s may be obtained in an equivalent manner to their 
quantum analogs (see [IB]), but these expression are not required for the discussion here. 
Thus, since the final state at t and s after the reabsorption of the photon is the same as 
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the initial one before its emission we may write for the relevant P 1 at the end of the whole 
process of emission and reabsorption [13] 



The coefficient C(t, s) denotes the mentioned evolution from the initial state P / (g°,t 0) So) to 
the final one P\q N , tjy, SAr|g°, to, sq) and is found as in [13] which discusses the same process 
in quantum terms (without using the variable s). We first note that the entire interaction 
of (emission+reabsorption) in the variables t and s, which is described after Eq f l24l ). may 
be written as a sum of two separate terms P(t) and P(s) each of them involves only one 
variable. These two terms are 



where we have set, as remarked, s = t = 0. Each of the two expressions P(t) and P(s) is, 
actually, an account of the whole process of emision and reabsorption, as discussed after Eq 
( 1241) . in the respective variables t and s. C(t, s) from Eq ( 1251) is given by the sum P(t) + P(s) 
so that in the equilibrium state obtained in the limit in which all the values of s are equated 
to each other the term P(s) vanishes and remains only the term P(t) as should be [TB]. This 
is because we have already equated the initial so to zero so for equating all the other s's 
to each other one have to set also the other values of s equal to zero which causes P(s) to 
vanish (see the second of Eqs (1261) ). Thus 



P I (q N , t Nl s N \q°, t 0l so) = P 1 {q° , t , s„) + C(t, s)P I (q°, t , s Q ) 



(25) 




(26) 




exp(i(e 2 + iS(e 2 + e 1 + w x ) +w x - e 1 )s 2 )ds 2 , 



C(t,s) = P(t) + P(s)=Y,9x s gi / dt 2 



(exp(z(ei - e 2 - w x )t 2 ) - 1) 



A. 



H e i - e 2 - w\ 



■ exp(z(e 2 + w x - ei)t 2 ) + ^ g Xs gt s ■ 



(27) 
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ds 2 r, —77 ; ; « exp(z(e 2 - ei + w x + 

. . \ \ \ \- gA^ r. / (exp(2(e 2 - ei + w A )t) - 1) 
x 8 H e i - e 2 - wa) H £ 2 - ei + io A ) 

+ y OXsOt ^2S(e 2+ e 1+W ,)s _ ^ 

x s ^( e i - e 2 - ^a s + i<5(e 2 + ei + wa)) «25(e 2 + ei + w A ) 
(exp(z(e 2 - ei + w x + iS(e 2 + e x + w A ))s) - T 



z(e 2 - ei + wa s + ^(62 + ei + w A )) 



-)] 



The first quotient in the square parentheses of the second sum, which is of the kind jj, may 
be evaluated, using L'hospital theorem [21], to obtain for it the result of s so that Eq f|27l) 
becomes 

C(t, -) = P(t) + P(s) = £ 7 ^ A t - { (eMi(e*-e 1+ wx)t)-l) 

x s < £ i - e 2 - ™a) «(e 2 - ei + w A ) 

+ Y" g^jk* [ s _ (28) 

^ i(ei - e 2 - w As + z5(e 2 + d + itf A )) 

(exp(i(e 2 - ei + w x + i5(e 2 + e 1 + wa))s) - 1) • 



i(e 2 -e x + w Xs + iS(e 2 + e 1 + w x )) 



-)] 



The last expression for C(t,s) yields terms of several kinds among them those which are 
proportional to t and s, others which are oscillatory in these variables, and also constant 
terms. Thus, for large t and s the oscillatory as well as the constant terms may be neglected 
compared to t and s as in the analogous quantum discussion of the same process [13]. 
Substituting the resulting expression in Eq (f25j) one obtains 



P I (q N , t N , s N \q°, 0, 0) = PV, 0, 0)(1 + C(t, s)) = P 7 (g°, 0, 0)(1 + itAe x + isAe Xa ), (29) 
where, 

Ae A = V 9Xs9 ^ , Ae As =V ^7 V (30) 

^e 2 -ei + w x x s e 2 ~~ e i + w *s ~ *<H e 2 + ei + w x ) 

The result in Eq (|29|) is only for the first-order term in Eq (l24j) . If all the higher order terms 

of this process are taken into account one obtains, analogously to the quantum analog (in 
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which the variable s is absent), the result 

P I (q N , t N , s N \q°, 0, 0) = PV, 0, 0)(1 + C(t, s)) = P 7 (g°, 0, 0)(1 + (itAe x + 

+ ^(ttAe x ) 2 + . . . + ^tAe x ) n + ■■■) + (isAe Xs + ^(itAe Xs ) 2 + ... (31) 

. . . + -^(itAe Xs ) n + ...) = P\q°, 0){e itAe > + e isAe ^ - 1) 

The right hand side of Figure 3 shows the diagram of the fourth order term of this process. 
Now, as required by the SQ theory, the stationary situations are obtained in the limit of 
eliminating the extra variable s which is done by equating all the s values to each other and 
taking to infinity. Thus, since, as remarked, we have equated the initial s$ to zero we must 
equate all the other s values to zero. That is, the stationary state is 

Urn P I {q N ,t N ,s N \q°, 0,0) = lim P 1 (q° , 0, 0){e itAe > + e isAe > - 1) = F J (g°, 0)e ltAe * (32) 

s— >0 s— >0 

The last result is the one obtained in quantum field theory [13] for the same interaction of 
(emission+reabsorption). The quantity Ae x , given by the first of Eqs (l30l) . has the same 
form also in the quantum version [EH fl4] . where it is termed the energy shift. This shift 
have been experimentally demonstrated in the quantum field theory for the case of a real 
many-state particle in the famous Lamb shift of the Hydrogen atom |T3l [HI [15]. 

As remarked, this electron-photon interaction is discussed in the literature [EH [HI [15] 
without using any extra variable and the result at the right hand side of Eq ([32l is obtained. 
This result have been obtained here by discussing the simulation process of the electron- 
photon interaction, that is, by using the extra variable s in the limit of equating all its 
values. In other words, as for the harmonic oscillator case, introducing the expression of 
the detailed electron-photon interaction into all the relevant codes (into the subintervals of 
s and t of all the programmers) yields a correlation among them which truly represents the 
corelation of the real interaction. Thus, as for the harmonic oscillator example in which the 
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use of the SQ method and the extra variable s lead to the real known quantum correlation 
(compare the two Equations (TT7l) . (fl8l ). so also here we obtain, using the same method and 
extra variable, the known expression for the probability P 1 (q N , t\q° , to) . Also, the writing 
stage of the relevant programs, which is characterized by the different versions of the code, 
is represented by the last result of Eq (13TI ) (the analogous writing stage of the harmonic 
oscillator example is given by Eq (JT6l)). The equilibrium stage corresponds to the case where 
all the values of s are equated to each other as seen from Eq (l32l) in which all the s's, 
including the initial so (see the discussion after Eqs (1241 ) and (I3T1)). are assigned the value 
of zero. The analogous expression for the harmonic oscillator example is given by Eq (fTTl) . 
The corresponding similarities between the correlations for the two remarked, 
the ultimate tests which show that the relevant simulations truly represent the simulated 
phenomena. 



Concluding Remarks 

We have discussed the code-writing stage of the numerical simulations of real phenomena 
using physical terms and terminology. The method applied for discussing this code-writing 
stage is the stochatic quantization method of Parisi-Wu-Namiki [7, 8J where an extra variable 
is introduced that takes account of an assumed stochastic process (in this variable) which 
allows a large number of possible different behaviours of the system. The equilibrium con- 
figuration is obtained [3, [8] when this variable is eliminated through equating all its different 
values to each other and taking to infinity. This equating of all the possible s values to each 
other introduces an element of repetitions of the same process through which the system is 
stabilized and brought to its equilibrium configuration. 

We discuss the system of a large ensemble of programmers which all try to simulate the 
same phenomena by writing the appropriate code which naturally will not be the same for all 
of them. The mentioned equating of all the s values means that all of them write the same 
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code and this constitutes a repetition of of writing the same program so that the probability 
to see the same numerical simulation is obviously large. 

We note that obtaining numerical equilibrium configuration through running a large 
number of times the related code upon the computer screen is the main characterstic of many 
simulations processes especially those concerned with finding numerical solutions of physical 
situations. For example, any one who try to numerically solve any differential equation 
which governs the evolution of some physical system knows, as shown in the following, that 
the solution suggested by the computer is obtained only after repeatedly updating the given 
differential equation. Better statistics is obtained when the number of iterations grows since 
this increases also the number of samples. The advantage of these repetitions is clearly 
seen for the case of simulating the long range correlation functions [221 [23] for which the 
conventional Monte Carlo simulation methods to numerically simulate them, using path 
integrals, fails [23]. It has been shown explicitly by Parisi [22], and Namiki et al [23], using 
SQ methods, that the following two point connected correlation function [16] 



d . [ d\q]qie h . . 1 . , . . . 

<M> = 55 ( 7lhr)l« = s«*>» - <«» + «(*>. P3) 

where the action Sh = S — hqi involves a small external additional term — hqi, may be solved 
by replacing the ensemble averages with and without the small external source < q\ >h and 
< qi > respectively by the time averages calculated from the following Langevin equations 

9S 

qi = + m 34 

oqi 
9Sh 

Qi = — + Vi (35) 
oqi 

The T) and fji are independent and assumed to satisfy 



<r] >=<r] l >= 0, <Vi(t)v m ( s )>=<Vi(t)Vm( s )>= 25i m 5(t- s) 
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Solving the right hand side of Eq ([331) generally results in a large statistical error [22, 23j, 
so Parisi [22] uses the same random forces in Eqs (1341) . (l35l) . that is, rj = fj which reduces 
considerably the statistical error as shown in [23] (see the Appendix there). Thus, in order to 
be able to simulate and obtain the long range correlation functions one must equate to each 
other not only all the different values of s which are related to r] and fj but also to equate rj to 
fj so that in the stationary limit one obtains, as remarked, the sought-for simulations. But as 
noted by Namiki et al [23] the last method, although works well for the fixed potential and 
the o(4) model [16, 20J, breaks down when one uses it to obtain the long range correlation 
function for the o(3) model possibly due to its large degree of nonlinearity. This situation 
is avoided in [23] by initiating a new round of repetitions where each one begins from the 
final configuration of the former. That is, in order to improve the statistical results one have 
first, as remarked, to increase the number of samples which is obtained by parallel updating 
of Eqs IJML f l35ll without and with the external source respectively using the same random 
forces for rj and fj. These steps which are sufficient, as remarked, for the fixed potential and 
the o(4) models end in a breakdown of the simulation for the o(3) model when the updating 
process continues. Thus, one must [23] stop this updating before break-down occurs and 
restart the whole procedure from switching again the external source and updating Eqs 
( 1341) , ( 1351) starting from the last stopped configuration as the initial one of the new round 
of updating. In other words, by only repeating the switching and the updating process one 
obtains, numerically, the sought-for long range correlation functions for the o(3) model. 

Moreover, it has been shown (24j [25], [26] , that these repetitions not only lead to numer- 
ical stabilization but when they are really performed (not just through clicking upon the 
computer keyboard) lead to a real stabilization. This phenomenon, termed the Zeno effect, 
have been validated both theoretically [24], [26], and experimentally [25J. The main charac- 
teristic of this effect is the preservation, through a large number of repetitions of the same 
measurement, of an initial state of the sytem (24] [25], or guiding it through a prescribed path 
of evolution, from a large number of possible paths [26] . 
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The principle of repetition have been shown in Section 3 for the Harmonic oscillator 
example (see Eqs (fT6l) -(rT7l) and the discussion there) where we see that when the same 
version of program is shared among the ensemble members then the probability to find the 
same simulation of the harmonic oscillator in all the screens is large. Moreover, when the 
number of times of performing this simulation, which is related to the number of subintervals 
(see the discussion after Eq (EES))) of the finite total paths in s and t of each member of the 
ensemble, becomes large so that the duration of each is small then the remarked probability 
is unity (see Eqs (fT9l) - (j20l) ). This is so since each member of the ensemble have in this case 
exactly the same code related to the Harmonic oscillator and the correlation among them 
is, therefore, maximal. The same state of affairs have been found also for the energy shift 
example discussed in Section 4. In this case the required correlation is obtained through 
summing the relevant Feynman diagram to all orders. This influence of repeating the same 
experiment a large number of times have been shown to be effective also for classical systems 
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Figure 1: The harmonic oscillator correlation function from Eq (TT71) as a function of the time 
t for the values of m = 1 and initial eigenvalue of ifo = 0.4. It begins from an initial value 
of 1.25 (which corresponds to the remarked values of m and w ), proceeds to a maximum 
value from which it descends to zero for large t. 




Figure 2: The harmonic oscillator correlation function from Eq (TT6l) as a function of the time 
t and the variable s for the same values of m = 1 and w = 0.4 as in Figure 1. The integral 
in Eq (l20l) have been numerically calculated for values of t and s in the ranges 1 < t < 20 
and 1 < s < 20. Note that the correlation tends to zero for large s even at those values of t 
in which the correlation from Eq (12TI ) (without s) obtains its larger values. 
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Figure 3: The left hand side of the figure shows the process of emitting and reabsorbing 
a photon in the time interval (to,t) where the energy is not conserved. The electron is 
represented in the figure by the directed arrow and the photon by the wavy line. The right 
hand side of the figure shows the same process repeated four times, in a perturbative manner, 
over the same time interval. 



